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Background: Spurious zero-energy Nambu-Goldstone (NG) modes appear when the symmetry of a system is 
spontaneously broken. The Thouless-Valatin inertia, the collective inertia of the NG mode, contains important 
information concerning collective motion. 

Purpose: To establish an efficient and precise method for deriving the collective inertia and the conjugate 
operator of a NG mode, we derive an expression for the response fnnction in terms of the coordinate-momentum 
representation of the quasiparticle random-phase approximation which is valid even if a symmetry-restoring zero- 
energy mode is present. 

Methods: We use the finite amplitude method for computing the response function of superfluid nuclei with the 
nuclear density functional theory. 

Results: We derived analytically the collective inertia and the conjugate coordinate operator of the NG mode 
from the zero-energy linear response with the momentum operator of the NG mode. The formulation is tested 
in the cases of translational and pairing rotational modes. Illustrative calculations are performed for the neutron 
pairing rotation in Sn isotopes, the proton pairing rotation in V = 82 isotones, and the neutron and proton 
pairing rotations around the ^®°Xe nucleus. 

Conclusions: The proposed formulation allows us to compute the collective inertia of the NG mode precisely and 
efficiently. The conjugate coordinate operator can be utilized to remove spurious contributions to the strength 
distribution in the finite amplitude method. 

PACS numbers: 21.60.Jz, 21.60.Ev, 21.10.Dr 

I. INTRODUCTION 

Spontaneous symmetry breaking is a universal phe¬ 
nomenon that plays an essential role in various fields of 
physics. The emergence of pions from chiral symmetry 
breaking and gauge-symmetry breaking in superconduct¬ 
ing systems are typical examples. 

The nucleus is a finite quantum system, whose ex¬ 
act ground state does not break the symmetries of the 
Hamiltonian. However, if we introduce a onc-body mean- 
field approximation such as the Hartree-Fock-Bogoliubov 
(HFB) approximation to the low-energy nuclear many- 
body problem, spontaneous symmetry breaking can take 
place to account for more correlations within the one- 
body approximation [1-5]. Continuous symmetries that 
are conserved in the nuclear Hamiltonian which can be 
broken spontaneously are the translational, rotational, 
and particle-number gauge symmetries. The isospin sym¬ 
metry can be broken spontaneously [6, 7], but is explicitly 
broken in the level of the effective nuclear and Coulomb 
interactions [8]. 

If a continuous symmetry is broken spontaneously, a 
zero-energy Nambu-Goldstone (NG) mode appears which 
restores it [9, 10]. In the case of the symmetry-broken 
nuclear mean field, the NG mode appears in the self- 
consistent solution to the quasiparticle random-phase ap¬ 
proximation (QRPA) as a symmetry restoration mode 
[11]. The NG modes that correspond to broken sym¬ 
metries in the mean field approximation (for example: 
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the center of mass mode for translational symmetry, the 
rotational mode for rotational symmetry, and the pair¬ 
ing rotational mode for particle-number gauge symme¬ 
try) are associated with infinitesimal transformations of 
the frame of reference. As such, NG modes are known as 
“spurious modes,” since they do not represent a physical 
excitation within the intrinsic frame. 

The actual spectroscopy measurement is performed in 
the laboratory frame for a nuclear energy eigenstate that 
preserves the symmetries of the Hamiltonian. Then the 
inertia of the NG mode in the intrinsic frame has an 
experimental correspondence. The collective inertia for 
the NG mode within the QRPA framework is called the 
Thouless-Valatin inertia [12]. The meaning of this in¬ 
ertia changes depending on the NG mode present: for 
translational motion it represents the total mass of the 
nucleus; for the rotational and pairing-rotational modes, 
it represents the rotational and pairing-rotational mo¬ 
ments of inertia, respectively. Except for the center of 
mass mode where both the coordinate and momentum 
QRPA phonon operators are known [11], the Thouless- 
Valatin inertia is not known in advance, and the QRPA 
equations must be solved in order to find it. 

The Thouless-Valatin inertia is important because it 
contains information on symmetry restoration, namely, 
ground state correlations to a particular broken sym¬ 
metry [13, 14]. This inertia was recently used for the 
three-dimensional spacial rotational mode in the five¬ 
dimensional quadrupole collective Hamiltonian [15, 16]. 
When compared with the Belyaev moment of inertia 
[17, 18] (its simplihed version), the Thouless-Valatin in¬ 
ertia is typically 30% higher; this is due to the contri- 
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buttons from two-body residual interactions not found 
in the Belyaev moment of inertia [19]. The Thouless- 
Valatin inertia has also been applied to symmetry break¬ 
ing in quantum dots [20]. Given its usefulness, a micro¬ 
scopic derivation of the inertia based on the nuclear en¬ 
ergy density functional (EDF) would be helpful in mak¬ 
ing systematic calculations with predictive power [21]. 

The expression for the Thouless-Valatin inertia is well 
known in terms of the QRPA A and B matrices [11, 22- 
24]. However, its direct evaluation is not trivial due to 
the size of these matrices. Therefore, several alternative 
approaches for computing the Thouless-Valatin inertia 
have been proposed: the cranked mean-field calculation 
[19, 25, 26] and a perturbative expansion based on adia¬ 
batic time-dependent HFB [27]. 

An alternative approach to solving the QRPA equa¬ 
tions based on linear response theory has been developed 
for nuclear density functional theory: it is called the finite 
amplitude method (FAM) [28]. In the FAM, the response 
function to an external one-body field can be evaluated 
only through a one-body induced field, and the calcula¬ 
tion cost is significantly reduced compared with the ma¬ 
trix diagonalization of the QRPA equations. The FAM 
has been implemented on various versions of Skyrme- 
HFB codes [29-31] and within covariant density func¬ 
tional theory [32, 33]. Not only does the FAM serve as 
an efficient way to calculate strength functions [34, 35], 
it has also been used to improve upon solution methods 
of the QRPA equations [36-38]. Up to now the FAM has 
been used to find physical excitations, and has not been 
formulated for the symmetry restoration NG modes. 

The NG mode sometimes gives an unphysical contribu¬ 
tion to a physical excitation that has the same quantum 
numbers. Its energy can be non-zero when a self consis¬ 
tency between the HFB and QRPA equations is broken, 
or when basis truncation violates the exact symmetry nu¬ 
merically. In Ref. [28], a prescription to remove such con¬ 
tamination from the FAM response function is proposed 
and demonstrated in the case of center of mass motion. A 
similar prescription is employed in the iterative Arnoldi 
diagonalization [39]. To apply this prescription to other 
modes, however, we need full solutions of the NG mode, 
which are not known in advance other than in the case 
of center of mass mode. 

The aim of this paper is to present a formalism of the 
FAM that can be applicable in the presence of the NG 
mode and give expressions for the Thouless-Valatin iner¬ 
tia and conjugate coordinate operator of the NG modes. 
The expression of the FAM response function in terms of 
the QRPA solutions [37] is based on the XY representa¬ 
tion of the QRPA equation, and is not applicable to the 
NG mode. In this paper, we derive the spectral repre¬ 
sentation of the response function of the FAM in terms 
of the momentum-coordinate (PQ) representation of the 
QRPA [11, 40]. This is done to handle both the NG 
modes and imaginary solutions of the QRPA in addition 
to the physical modes. 

This paper is organized as follows. In Sec. H, the 


QRPA equations in the PQ representation are recapit¬ 
ulated. Section HI contains a brief introduction to the 
FAM, and its response function is expressed in the PQ 
representation of the QRPA in Sec. IV. Then in Sec. V, 
the expression for the Thouless-Valatin inertia is derived. 
Numerical tests of the formalism for the center of mass 
motion and pairing rotations are presented in Sec. VI, 
and realistic calculations of the neutron and proton pair¬ 
ing rotations are shown in Sec. VH. Lastly, conclusions 
and outlook are given in Sec. VHI. 

II. QRPA IN XY AND PQ REPRESENTATIONS 

We recapitulate the QRPA matrix expression for later 
convenience [11]. The QRPA equation in XY represen¬ 
tation is given by 

.SA’ = T. 3 XO, (1) 


where 



where are the two-quasiparticle indices, ij are 

the indices for QRPA eigenmodes, A and B are the 
QRPA matrices, and Hi is the energy of the QRPA eigen¬ 
mode. The indices of E3 can be either fj,iy, or ij. The 
matrix Q in O is flij = fliSij. The amplitudes A* and V* 
are the two-quasiparticle amplitudes of a QRPA phonon 
operator 

ol = E (4) 

with they are normalized with 

{[4,6]])=^.,, {[6J,6]]) = o, (5) 

or equivalently 

AtEgA = S 3 , ASgAt = Eg. (6) 

The QRPA equations are written in the Hamiltonian ex¬ 
pression as 

[-^QRPA, 4] = [HQKPA^Oi] = (7) 

where Hqbj>a is the QRPA part of the Hamiltonian. We 
always have pairs of solutions with positive and nega¬ 
tive frequencies. This expression assures that the two- 
quasiparticle part of the Hermitian broken-symmetry op¬ 
erator Png given by 

Png = (8) 

lJs<y 
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is always a solution Ong of the equation with IIng = 0 - 
The problem with the XY representation of the QRPA 
equations for NG modes is that they cannot be normal¬ 
ized with Eq. ( 6 ). 

We now switch to the PQ representation to handle NG 
modes. We introduce Hermitian coordinate and momen¬ 
tum operators that describe the eigenmodes as 

4 =V^^(4 + 0 \) = ^ 

(9) 

V, _ 6t) = ^ 

fl<l' 

( 10 ) 

and regard these operators as a starting point. Here Mi 
is the inertia for each mode. The QRPA equations in 
terms of these operators are written as 


where Mng is the Thouless-Valatin inertia for the NG 
mode. In the PQ representation, Png and Qng are nor¬ 
malizable with Eq. (15) and the Thouless-Valatin inertia 
is given by 

2Png(A -I- P)“^Png (IhiPng = 0), 
-2Png(A- P)-iPng (ReP ng = 0 ). 

( 20 ) 

The expression above depends on whether Png is real or 
pure imaginary. Computation of the Thouless-Valatin in¬ 
ertia from Eq. (20) requires full evaluation of the A and 
B matrices of large dimensions for recent nuclear den¬ 
sity functional theory without symmetry restrictions. In 
the next section, we derive expressions for the Thouless- 
Valatin inertia based on linear response theory. 

III. FINITE-AMPLITUDE METHOD 



E 




M. (^0,/ 

( 12 ) 


In this section we introduce the FAM, and express the 
response function in the PQ representation following the 
notations in Ref. [30]. We start with an external time- 
dependent field F{t) with a frequency uj and a small finite 
amplitude parameter rj 

P(t) =r;|Pe-“‘-f Pte“‘| (21) 


By defining the matrices 


V, 




P^ 

_pi* 


Q 

Q 




( 0 

0 

(13) 


applied to the system. Here the one-body operator is 
written in the quasiparticle basis 

^ = E ^ P“S^., (22) 

/L<L' 


with Mij = MiSij, the QRPA equations in PQ represen¬ 
tation are summarized as 

SV = EaVW. (14) 

The operators Pi and Qj are normalized with 


where In the FAM, we solve the time- 

dependent HFB (TDHFB) equations with the external 


field 





(23) 


([Si, Fj\) —iSij, 

or equivalently 


([Si, Qj]) = il'Pii'Pj]) =0, (15) dependence is governed by the forced oscilla¬ 

tion of the external field. The time-dependence of the 
quasiparticle and the Hamiltonian is given by 


V^Y3V = Y2, VS2Vt = E3, 

where 



(16) 


(17) 


The solution for a zero-energy NG mode is written as 


[A PW Png \ 
[b* a*) {-P*J 

(A B\ f Qng \ 
[b* a*) [-q*^J 


= 0 , 


i ( Png\ 

M-^g V NG/ ’ 


(18) 

(19) 



(24) 

Sa^it) =vJ2^t {V,^(o.)e—* +y;4cc)e“‘} , 

1/ 

(25) 

H{t)=Y^E^B^^ + 5H{t), 

(26) 

SH{t) =77 {<5P(a;)e-*‘^‘ -f <5Pl'(a;)e“‘} , 

(27) 

SH{u;) = + SH^^Uoj)A^.] , 

(28) 

where is the quasiparticle energy and and 

Y^y{uj) are the FAM amplitudes. Using the expressions 
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above and taking the terms linear to rj (small-amplitude IV. FAM RESPONSE FUNCTION IN THE PQ 
approximation), the TDHFB equations are written as REPRESENTATION 


{E^ + E,- oj)X^,{u;) + (29) 

(E^ + E,+ w) V(c^) + = - F°l. (30) 

This expression does not involve the A and B matri¬ 
ces explicitly, but through the one-body induced fields, 
6E[^^{uj) and SE[^^{u}). This is the advantage of the FAM: 
the response of the system to the external field F can be 
evaluated from one-body quantities only. We refer more 
details on how to compute the one-body induced field in 
nuclear density functional theory to Ref. [30]. We also 
note that the FAM equation is formally equivalent to 
linear response theory, and can be written as 



[A B\_ (I 0\ 
\B* A*) ^1^0 -l) 



where R{uj) is the response function. 



(31) 


Now we derive the expression of the response function 
R{u>) in terms of the PQ representation of the QRPA. A 
similar derivation has been done for the XY representa¬ 
tion of the QRPA in Ref. [37]. Using Eqs. (14) and (16), 
the response function i?(w) is written as 

R{uj) =[S - wEa]-! = V[W - a;J]-^E2V^ (32) 


where 




( —uj —iM 

-UJ ) ’ 


(33) 


and I is a unit matrix. We note that this expression 
cannot be defined at w = This excludes w = 0 if a 
NG mode is present. This matrix has blocked structure 
for the same index in the four blocks, and we can take 
the inverse for each 2x2 matrix. The response function 
and the FAM amplitudes with P and Q coefficients are 
now given as 


R{u:) 


E' 




(lujPl, - )QU,, + (-M-^P» - lUjQUPX, )QU,, + + ^UJQ 


{-lujPll + M.nlQi 


.)Q 


li' u' 


+ {M-^Pll+iujQ^;,)Pl1,, 


X, 


t,u{uj) 

(w) 


Y 


E/ 

fi'v 


^1/ I fl'u' 

{lujPll - M,n^Q»)QY, + i-M- 


^Pll - lUjQ 


(uj) 


p20 
a'v' 
f02 
y'. 


E 


1 


LO 




(-zccP;, + M,u2Q;j(g,iPio) 
(luupn - )(Q,1F10) + 


—p* 
1 


iujQ\ 





where we define the transition strengths from the ground strength function is 
state to the states expressed with Vi and Qi phonon op¬ 
erators as 


(PIPJO) ={[p, P]) = ^ P%Ff, - P;,F°l, (36) 

^<y 

(g.lPjo) =([ 4 , F]) = Q%Fll - QluFll- (37) 


^<y 


S(F,uj) = Y FlTX^,(uj) + Fll*Y^,(uj) 

fi<y 

= E + M.n^m\F\0)\^ 

i I K 

+a;[gP],(P)}, (38) 


We note that the FAM amplitudes X{uj) and Y(uj) are 
not A* and U* eigenvectors of the QRPA themselves, 
therefore they are well-defined through the linear re¬ 
sponse equation even if there are NG modes. The FAM 


where we define a real quantity \QP]i(F) = 
z((g,lPl0)*(P|P|0)-(PlPl0)*(g,|P|0)). when NG 
modes are not present (fli ^ 0), the following transition 
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strength can be defined: 

and by substituting these into Eq. (38), we can go back 
to the original expression of the FAM strength function 
[37] 


to the strength distribution. The procedure to remove 
this contribution from the center of mass mode has been 
proposed in Ref. [28]. 

V. THOULESS-VALATIN INERTIA FOR 
NAMBU-GOLDSTONE MODES 

A. Thouless-Valatin inertia from the momentum 
operator 




( l(^l^| 0 )P 
^ V ^ 


Ko^\ 

UJ I 


(41) 


Generally, the solutions of the QRPA equations consist 
of physical modes with Of > 0, NG modes with Of = 0, 
and imaginary modes with Of < 0. Imaginary solutions 
of the QRPA equations can occur when the HFB state 
does not correspond to a variational minimum. The PQ 
representation of the strength function (38) is valid for all 
three kinds of the modes. We consider the QRPA equa¬ 
tions at a non-variational minimum mainly in two cases. 
When the HFB code has a symmetry restriction, we can¬ 
not take the variation against the restricted degrees of 
freedom, and unexpectedly the HFB state obtained with 
the symmetry-restricted code can be unstable. Typi¬ 
cal examples are when a deformed nucleus is computed 
with a spherical HFB code, or when a triaxial state is 
computed with an axial HFB code. The transition to 
isoscalar pairing condenstation with a HFB code with 
proton-neutron symmetry and proton-neutron particle- 
particle RPA was recently discussed from this point of 
view [41]. Another case is large-amplitude collective mo¬ 
tion. The local QRPA [16] based on the adiabatic theory 
of large-amplitude collective motion [42] requires the so¬ 
lutions of the QRPA equations at non-equilibrium HFB 
states. 

We can split the contributions of three modes to the 
strength function as 

SiF, Lu) =SiF, w)phys + SiF, w)ng + S{F, w)i„,ag, (42) 


where the contribution from the NG mode is given by 


The momentum operator T^ng of a NG mode is a con¬ 
sequence of a broken symmetry of the system. By using 
it as an external field of the FAM, we can show that the 
contribution to the strength function is zero 

^(Png,w) = 0 (45) 


from Eqs. (15), (38), and (43). Here we recall that the 
response function (34) is undefined at w = 0 in the pres¬ 
ence of the NG mode. From Eq. (31), the linear response 
equation at w = 0 is written as 


(A B\ (X{0)\ 

[b* a*) [¥{0)) [P^J ■ 


(46) 


Assuming that the A and B matrices are real, we have 


(A + H)[A(0) + r(0)] =-2 Png (ImPNG = 0), (47) 
(A-P)[A(0)-T(0)] =-2Png (RePNG = 0). (48) 
The FAM strength function at a; = 0 is then 


5(Png,w = 0) = ^(Pi;g)m^^m-^(0) + {Png)>^M0) 

_ f —2Png(^ +-S)“^Png (IniPNG = 0) 
'|2Png(^ —-S)“^Png (R.ePNG = 0) 

= — -Mng, ( 49 ) 


where we use Eq. (20). Therefore the FAM strength func¬ 
tion for the momentum operator of the NG mode is sum¬ 
marized as 


5'(Png, w) 


0 iy ^ ifli) 

—Mng (w = 0 ) 


(50) 


s'(p, w)ng = 

i,ni=o 


MiUJ^ 


[QPUn 1 


(43) 


and the strength at zero frequency gives the Thouless- 
Valatin inertia. The coordinate opeartor of the NG mode 
is then given by 


The same expression is found in Ref. [43]. When we com¬ 
pute the strength function distribution with the FAM, we 
replace the real frequency w by a complex value for the 
frequency uj+ij where the imaginary part gives the width 
r = 27 . If the external field can excite the NG mode, 
there is a spurious contribution 

--ImS'(P,w -I- i7)NG 

TT 

^7 ^ f 2c^|(P.|P|0)p [QP].(P) 1 

a;2 + 72 J ^ ^ 


Qng 


*-^Ng(^ + P) ^^NG (IhiPnG = 0 ) 
-iM^Q{A — P)“^Png (Re Png = 0 ) 


r A(o) + y(o) 

I 25(Png,0) 
1 . x(o)-y(o) 

[ 25(Png,0) 


(ImPNG = 0), 
(Re Png = 0). 


(51) 


Except for the trivial case of the center of mass mode, 
the coordinate operator of the NG mode is not known in 
advance, and this expression will be useful for removing 
spurious modes [28]. 
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B. Thouless-Valatin inertia from the coordinate 
operator 

An alternative derivation of the Thouless-Valatin iner¬ 
tia is found from the FAM calculation using the conjugate 
coordinate operator of the NG mode Qng as an external 
field. In this case, the strength function is derived from 
Eqs. (15) and (43) as 

S{Qng,uj) = (w ^ ±^^i)- (52) 

Therefore the Thouless-Valatin inertia is given from the 
energy-weighted sum rule [11] as 

Mj^q = 2mi(QNG) = [ wS'(Qng, w)dw, (53) 

27’"* JAi 

where Ai is the counterclockwise half circle in the com¬ 
plex energy plane from uj = —IRai to iRai centered at 
the origin [38]. 


C. Thouless-Valatin inertia and the m_i snm rnle 

The expression for the Thouless-Valatin inertia (20) 
is same as the inverse-energy weighted sum rule given 
in Ref. [44] but for the momentum operator of the NG 
mode. Moreover the dielectric theorem [11, 44] connects 
the inverse-energy weighted sum rule through the con¬ 
strained HFB state 1(/>(A)) computed with the Hamilto¬ 
nian H + A'Png 


7-i(^ng) = 


dX 


(<^(A)|75ng|</>(A)) 


A=0 


A4ng 

2 ■ 

(54) 


This theorem allows us to compute the Thouless-Valatin 
rotational moment of inertia from the cranked HFB cal¬ 
culation [19, 20, 25, 45]. These facts show that the 
Thouless-Valatin inertia is related to the inverse-energy- 
weighted sum rule m_i('PNG)- However we need a careful 
consideration of the m_i sum rule when NG modes are 
present. The expression of the m_i sum rule in terms of 
the transition matrix elements of an operator F is given 

by 


m_l(A) = ^H-ll(^lA|0)l^ (55) 

i 

and the NG modes should not be included in the sum¬ 
mation. Because of the discontinuous character of the 
FAM strength function S'('Png,w = 0), the inverse 
energy-weighted sum rule from the contour integration 
of Ref. [38] does not provide the Thouless-Valatin iner¬ 
tia. The contribution from the NG mode to the m_i sum 
rule has to be added separately. A similar discussion is 
found for the mi sum rule in the Appendix of Ref. [46]. 


D. Approximate symmetry 


Even though the Hamiltonian preserves its symmetry, 
it can be explicitly broken in the settings of a numeri¬ 
cal calculation. This is actually the case for translational 
and rotational symmetries. We often express the single¬ 
particle states either in the harmonic oscillator basis ex¬ 
panded about the center of mass or in the coordinate 
lattice of a finite box. Ideally the translational displace¬ 
ment leaves the energy of the system invariant. But if 
the wave function is expanded in a finite basis, the trans¬ 
lational/rotational symmetry is broken explicitly. This 
results in the translational and rotational modes appear¬ 
ing at finite energies. We discuss the effect on such ap¬ 
proximate symmetries on the evaluation of the Thouless- 
Valatin inertia. When the symmetry is approximate, the 
excitation energy Hng is finite, and the contribution to 
the strength function is actually the same as other phys¬ 
ical modes. The momentum and coordinate operators 
of the NG mode are now only approximate solutions to 
the QRPA equations. Therefore we have the following 
approximate expressions for the strength function at en¬ 
ergy around w = Hng: 

(56) 

^(iPNG,u;)- ^^°^r . (57) 

At w = 0, the contribution from the finite excitation 
energy Hng in Eq. (57) is canceled, and Eq. (49) remains 
as a good approximation to the Thouless-Valatin inertia. 
The energy-weighted sum rule (53) is valid if Rai > AIng- 
The position of the excitation energy of the NG mode is 
estimated from Eqs. (56) and (57) as 


n 


2 _ 
NG — 


1 

SiVyiG, 0)S'(QnG, 0) 


(58) 


E. Inglis-Belyaev cranking inertia 


The Inglis-Belyaev cranking inertia [17, 18] is easily 
evaluated when the FAM routine is available. The gen¬ 
eral expression of the cranking inertia for the Hermitian 
operators Fi and Fj is given by 

T7'i20* 77'720 _i_ 77'220 Z7'720* 

M,e(.a = 2y: ^-- . (ss) 


By setting the induced field and the energy to zero in 
the FAM calculation = 0 and w = 0), the 

FAM amplitudes for operator F) is given without self- 
consistent iteration: 


= 0 ) 

Yf,^{F,,u} = 0 ) 


pi ,20 

EfJ, + Ey 

pz,02 

En + El, 


(60) 

(61) 












The Inglis-Belyaev inertia is given from the amplitudes 
for another operator Fj as 
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(62) 


If Fi = Fj, the Inglis-Belyaev inertia gives the response 
function without an induced field. This derivation agrees 
with the fact that the cranking inertia does not include 
the contribution from the residual interaction. 


VI. NUMERICAL TESTS 

To check the derivation in the previous section and 
to show its applicability, we performed numerical calcu¬ 
lations for ^®Mg with the EDF SLy4 [47] with volume 
pairing Vq = —125.2 MeV fm^ and a 60 MeV quasi¬ 
particle energy cutoff. The HFB and FAM calculations 
were based on the computer code HFBTHO [48, 49] and 
its FAM extension to the K = 0 mode [29]. The cal¬ 
culations were performed with various sizes of the har¬ 
monic oscillator space from Ngh = 5 to 20. We used 
the proton-superconducting oblate deformed HFB state 
with /? = —0.18 and Ap = 0.619 MeV. The broken sym¬ 
metries in this HFB state are translational, rotational, 
and proton particle number. 


A. Center of mass mode 

The center of mass mode appears as a translational 
symmetry restoration mode. This is the only case where 
both the coordinate and momentum operators are known 

1 ^ . 
i=l i—1 

and the Thouless-Valatin inertia is nothing but the total 
mass Mcm = mA [11]. We use the z-component [K = 0) 
of the center-of-mass coordinate operator (Qcm)z as an 
external field of the FAM. Figure 1 shows the FAM 
strength function with various sizes of the harmonic oscil¬ 
lator single-particle model space Agh- In this calculation 
a real value of w is used, therefore the strength function 
'S'((<5 cm)2 , is also real. As we discussed in Sec. VD, 
translational symmetry is not an exact symmetry of the 
system if the single-particle states are expressed in a fi¬ 
nite harmonic oscillator basis, and the strength function 
is approximated by Eq. (56). Therefore the value of the 
strength function at w = 0 is finite, and there is a low- 
energy pole that corresponds to a spurious excitation of 
the center of mass mode. Its energy approaches zero as 
the model space size increases. The strength function 
for Ash = 20 is close to that for Ash = 15, indicating 
that further convergence of the spurious energy to zero 
is numerically difficult. Figure 2 shows the same plot 
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FIG. 1. (Color online) The FAM strength function 
S'((Qcm)z, ta) for a response to the center of mass coordinate 
operator with a real frequency uj. 



FIG. 2. (Color online) The FAM strength function 
S'((Rcm)z,uj) for a response to the center of mass momen¬ 
tum operator with a real frequency uj. 


but for the center-of-mass momentum operator (Pcm)z- 
This curve is approximated by Eq. (57), and the value at 
u! = 0 shows the Thouless-Valatin inertia. 

In Table I, we list the Thouless-Valatin inertia for the 
center of mass motion in the form of I/2m = A/2Mcm- 
The value for SLy4 is 20.73553 MeV fm^. The one-body 
center-of-mass correction to the kinetic energy [50-52] is 
not taken into account in this calculation because it ef¬ 
fectively scales with the nucleon mass m, and makes the 
comparison complicated. In a smaller harmonic oscillator 
model space such as Ash = 5, the center of mass excita¬ 
tion has higher excitation energy. This indicates the co¬ 
ordinate and momentum operators expressed in a small 
basis are not a good approximation to the NG solution 
of the QRPA, and the component of the center of mass 
can be distributed over the physical modes. The normal¬ 
ization of the coordinate and momentum operators are 
satisfied within 0.2% accuracy. With a larger harmonic 
oscillator model space, the spurious energy of the center 















of mass reaches zero, and the normalization of the oper¬ 
ators is more accurate. In Table I the Thouless-Valatin 
inertia computed from the energy-weighted sum rule of 
the center of mass coordinate operator with Eq. (53) up 
to 2 MeV and from the strength function with the center 
of mass momentum operator at zero energy with Eq. (49) 
are listed. Both inertias are close to the exact value, 
but the inertia computed from the momentum operator 
agrees precisely, beacuse it is evaluated at a single point 
w = 0 in the complex-energy plane, and is more precise 
than the discretized contour integration with the coor¬ 
dinate operator. Practically, in all cases of symmetry¬ 
breaking modes other than the center of mass mode, we 
know only the momentum operator, and the Thouless- 
Valatin inertia can be derived from it. This analysis using 
the center of mass mode shows that even if the symmetry 
is approximate, the Thouless-Valatin inertia can be com¬ 
puted very accurately from the EAM strength function. 

In the same table, we also show the Inglis-Belyaev in¬ 
ertia for the center of mass motion. As is well known, 
the Inglis-Belyaev value of the inertia deviates from the 
Thouless-Valatin inertia, showing that the contribution 
of the residual interaction is very important even in the 
case of trivial center of mass motion. 


B. Pairing rotational mode 



FIG. 3. (Color online) The FAM strength function for the 
proton pairing-rotational angle operator S(0p,(jj) for a real 
frequency uj for ^®Mg, computed with a Ash = 5 model space, 
plotted together with the curve in Eq. (52) with the Thouless- 
Valatin inertia Mtv =1.1545 MeV“^. 

5 '((Qcm)2,w) in Fig. 1. In the case of the pairing rota¬ 
tional mode, the symmetry is exact, and the position of 
the pole is exactly at^zero energy. The Thouless-Valatin 
inertia from the mi(Op) sum rule through Eq. (53) with 
Rai = 1 MeV is 1.1545 MeV“^, and is perfectly consis¬ 
tent with the inertia from the momentum operator JVp. 


Unlike center of mass motion, particle number opera¬ 
tors are defined in the configurational space. Therefore 
the broken particle-number gauge symmetry is always ex¬ 
act if the same model space is employed in the HFB and 
QRPA frameworks. 

As we have seen in Eq. (50), the strength function 
obtained with the proton particle number held, that 
is a broken particle-number gauge symmetry in ^®Mg, 
S{Np,uj) = 0 is zero except at w = 0. The FAM enables 
us to compute the discontinous value at w = 0 without 
any convergence problems, and elsewhere the strength 
function is numerically zero. The Thouless-Valatin in¬ 
ertia for the proton pairing rotation in the present case 
with Ash = 5 is 1.1545 MeV“^. The same quantity but 
for the neutron number operator is exactly zero including 
at a; = 0, because the neutron particle number is a con¬ 
served symmetry, and commutes with all other QRPA 
modes of a particle-hole type. 

From Eq. (51), we can derive the proton gauge angle 
operator 0p, which is the conjugate coordinate operator 
of the proton pairing rotation. Although we can com¬ 
pute the Thouless-Valatin inertia from S{Np,u! = 0), 
to check the consistency, we show the response to the 
coordinate operator of the proton pairing rotation in 
Fig. 3. In the figure, the strength function computed 
from the FAM and Eq. (52) with the Thouless-Valatin 
inertia determined from the proton number operator are 
compared. Both curves agree very well as expected. We 
note that this strength function S{Qp,uj) looks very dif¬ 
ferent with that of the center of mass coordinate operator 


VII. REALISTIC EXAMPLES 
A. Pairing rotation in single-closed shell 

As for realistic examples of NG modes in nuclei, we dis¬ 
cuss pairing rotations in single-closed shell nuclei. The 
ground states of tin isotopes are known to form a neu¬ 
tron pairing rotational band [11, 53, 54], which has been 
actively discussed recently in connection with the two- 
neutron transfer reaction [55-58]. The binding energy of 
the A-neutron isotope B{N, Z) from the nearby reference 
state with neutron Aq is written as 

(AA)2 

R(A, Zo) = B{No, Zo) + A„(Ao, Zo)AN + 

(64) 

where A A = A — Aq, and A„(Ao,2'o) is the neutron 
chemical potential of the reference state. Under the as¬ 
sumption that the ground states of even-even nuclei form 
a pairing rotational band, the last term can be regarded 
as the pairing rotational energy, and the neutron pairing 
rotational moment of inertia Jn{No, Zq) at a reference 
state is computed with the FAM from the zero-energy 
response of the neutron particle-number field 

J-„(Ao,Zo) = -5(A„,w = 0) (65) 

at a reference state. 
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TABLE I. The inertia for the center of mass motion of in units of MeV fm^. The Thouless-Valatin inertia is computed 
from the energy weighted sum rule of the center of mass coordinate operator (Qcm)z with Eq. (53) with a radius Rai = 2 
MeV discretized with Nai = 12 points, and from the strength function at zero energy with the center of mass momentum 
operator (f’cM )2 with Eq. (49). The Inglis-Belyaev value of the inertia, the spurious excitation energy of the center of mass 
mode evaluated with Eq. (58), and the commutation of the coordinate and momentum operators are also listed. 


Wh 

l/2m from (Qcm)^ 

l/2m from (R:m )2 

Inglis-Belyaev 

flcM MeV 

([(Qcm) 2 , {PcM)z\)/i 

5 

20.69748 

20.74676 

26.04977 

1.346 

0.998836 

10 

20.78073 

20.82140 

25.87571 

0.889 

0.999310 

15 

20.73573 

20.73232 

25.73650 

0.151 

1.000026 

20 

20.73946 

20.73666 

25.74138 

0.146 

1.000041 

exact 

20.73553 

20.73553 

- 

0 

1 


Again we use SLy4 with volume pairing in a Nsh = 20 
model space. The pairing strength is fixed to Vq = 
— 178.81 MeV fm^ which reproduces the experimental av¬ 
eraged neutron pairing gap A^^^(^^®Sn) = [Ai^^(^^®Sn) -|- 
Ai^^(^^'^Sn)]/2 = 1.100 MeV [59]. This parameter setting 
is used for all the calculations in Sec. Vll. We calculated 
the neutron pairing rotational moment of inertia at the 
reference state ^^®Sn {N = 66) that is located in the 
middle of the N = 50 and 82 shell gaps. The Thouless- 
Valatin moment of inertia for the neutron pairing rota¬ 
tion was Jn,TV = 5.95 MeV“^, while the Belyaev inertia 
was given by Jn,m = 4.71 MeV“^, and in this case the 
residual interaction contributes to the inertia about 30%. 

Figure 4 shows the neutron pairing rotational energy 
measured from ^^®Sn as a function of neutron number. 
The Thouless-Valatin inertia for the neutron pairing ro¬ 
tation explains the pairing rotational spectrum especially 
well in the vicinity of the reference state. Agreement be¬ 
tween the pairing rotational energy from the Thouless- 
Valatin inertia and the HFB energy shows the validity 
of the dielectric theorem far from the reference state and 
at small anharmonicity. We note that inclusion of the 
one-body center-of-mass correction violates the dielectric 
theorem, because the correction term (1 — 1/A) is not 
variational with respect to the change of particle num¬ 
bers [60]. As the neutron number changes from N = 66, 
deviation from the parabola curve is seen, both in the 
HFB calculations and the experimental data, indicating 
a change of the intrinsic structures as a function of the 
neutron number. The deviation is larger in the neutron- 
deficient side. 

The next example is the proton pairing rotation in 
N = 82 isotones, where we take ^^^Nd {Z = 60) as 
our reference state. The averaged proton pairing gap 
A^3) (i42jv,jd) = 0.788 MeV is well reproduced with the 
same pairing functional parameters. We can define the 
proton pairing rotational moment of inertia Jp{Nq, Zq) 
by repeating the same discussion with Eq. (64) but for 
protons. The Thouless-Valatin inertia for this proton 
pairing rotation is 2.35 MeV“^, while the Belyaev in¬ 
ertia is 6.13 MeV“^. The residual interaction changes 
the pairing rotational moment of inertia by a factor of 
about 2.6. This is because of the residual Coulomb con¬ 
tribution: the Coulomb interaction is known to affect 



FIG. 4. (Color online) Neutron pairing rotational energy of 
even Sn isotopes. The red solid (blue dashed) curve is the 
pairing rotational energy {N — 66)^/[2 J^j^tvcb) (^^®Sn)] with 
the Thouless-Valatin (Belyaev) pairing rotational moment of 
inertia evaluated at ^^®Sn. The red squares are the HFB 
energy £;hfb(A)-£;HFB(“®S n) - A„.HFB(“®Sn)(A-66), and 
the black circles are the experimental values evaluated with 
— Bexp{N) — Xn,expC^^Sn){N — 66), where the binding energy 
is taken from Ref. [59]. The experimental neutron chemical 
potential is evaluated with A„,exp(^^®Sn) = [Bexp(^^®Sn) — 
Bexp(“^Sn)]/4 = -8.345 MeV. 


the proton pairing energy and gap at the mean-field level 
[61-63]. Because the Coulomb energy is proportional to 
Z^, its residual part directly contributes to the proton 
pairing rotational moment of inertia in the QRPA level. 
Figure 5 shows the proton pairing rotational energy mea¬ 
sured from ^^^Nd. To explain the experimental curve, the 
contribution of the residual interaction is essential in the 
case of the proton pairing rotation. 


B. Mixing of neutron and proton pairing rotational 
modes 

When both the neutron and proton are in the super¬ 
conducting phase, the neutron pairing and proton pair¬ 
ing rotational modes appear as NG modes. These zero- 
energy modes are degenerate because they have the same 
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FIG. 5. (Color online) Same as Fig. 4 but for the proton 
pairing rotational energy for even N = 82 isotones measured 
from {Z = 60). The experimental proton chemical 

potential is Ap,exp(^"^^Nd) = —5.761 MeV. 


K'^ = O'^' quantum numbers. Thus the eigenmodes as 
solutions of the QRPA equations are generally the linear 
combination of the two pairing rotational modes. The 
momentum operators of the two NG modes are written 
as 

Ni =Nn cos 9 + aNp sin 9, 

N 2 = — Nn sin 9 + aNp cos 9, (66) 

where 0 is a mixing angle, and a is a scaling parameter. 
The parameter a should be one under the isospin sym¬ 
metry, but we keep a general expression here. The overall 
scaling of the operators iVi and N 2 does not change the 
physics. The conjugate angle operators are written as 


where 


S{Nn, iV„) = - 2iV„(A + = S{Nr,,uj = 0), 

(70) 

SiNr,,Np)=-2Nr,iA + B)-^Np = SiNp,N^), (71) 
S{Np, Np) = - 2Np{A + B)-^Np = S{Np, uj = 0) (72) 


are the strength functions obtained from the FAM 
with zero frequencies. Here Nn and Np are the two- 
quasiparticle amplitudes of the particle number operators 
Nn and Np, respectively (we assume Im iV„ = Im iVp = 
0). The off-diagonal term S{Nn, Np) is obtained by eval¬ 
uating the strength function for the proton particle num¬ 
ber operator from the FAM equations with a neutron 
particle-number external field, or vice versa. The two 
parameters 9 and a are constrained from the orthogonal¬ 
ity of the two NG modes 


tan 29 = 


2aS{Nn,Np) 

SiNn,Nn)-a^SiNp,Np)' 


(73) 


In comparison with experimental data, we are interested 
in the pairing rotational energy in terms of neutrons and 
protons rather than in terms of the neutron-proton mixed 
eigenmodes. The pairing rotational energy is written as 


E,ot{N,Z) 


(AiVi)^ (AjV2)^ 

2 Afi 2 M 2 

A(AJY az)j- (“) 

(AA)2 ^ 2(AN)(AZ) ^ (AZ)2 
“ 2Jn„ 2Jnp ^ 2Jpp ’ 


(74) 


where ANi, AZi, AN and AZ are the deviation of the 
particle numbers from a reference state. The inertia ten¬ 
sor J is given by 


©1 =0nCOS0 -I-0pSin0 

a 

02 = — ©raSind-I-— 0pCOS0, (67) 

a 

where 0„ and Op are the neutron and proton two- 
quasiparticle parts of the operator. The operators in 
Eqs. (66) and (67) satisfy the commutation relations of 
Eq. (15) if the neutron and proton parts of the operators 
are normalized with the same commutation relation. 

The Thouless-Valatin inertias for the two NG modes 
are derived from Eq. (20): 

Ml = - S(Nn, Nn) cos^ 9 - a'^S{Np, Np) sin^ 9 

— 2aS{Nn,Np) sin 9 cos 9, (68) 

M 2 = - S{Nn, Nn) sin^ 9 - a^S{Np, Np) cos^ 9 

+ 2aS{Nn,Np) sin 9 cos 9, (69) 


-l ^_fSiNn,Nn) SiNn,Np)\ " ^ f 1 /Jnn 1 /Jnp\ 
\S{Np,Nn) S{Np,Np)J \1/Jpn 1/JppJ’ 

(75) 


and the neutron and proton components are explicitly 
written as 

^ S{Nn, Nn)S{Np, Np) - S{Nn, Np)^ 

\Jnn — , V 

S{Np,Np) 

„ S{Nn,Nn)S{Np,Np) - SjNn^Np)^ 

SiNn,Np) 

„ S{Nn, Nn)S{Np, Np) - S{Nn, Np)^ 

SiNn,Nn) 

This shows that the principal axes of the pairing rotations 
are not aligned with the neutron and proton directions 
in the gauge space because of the presence of the off- 
diagonal term S{Nn, Np). We note that this does not 
exist in the Belyaev inertia, because the two-quasiparticle 


(76) 

(77) 

(78) 
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indices in Eq. (59) are either neutrons or protons when 
neutron-proton mixing is absent in the mean field. The 
residual interaction plays an essential role for generating 
the neutron-proton term in the pairing rotational energy. 
In Ref. [64], the principal axes tilted against the neutron 
and proton gauge-angle space have been reported in the 
reduced energy kernel when neutrons and protons are 
superconducting. 


C. Neutron and proton pairing rotations around 

We consider the neutron and proton pairing rotations 
by taking a reference state at the open-shell deformed nu¬ 
cleus ^^°Xe. The lowest energy HFB solution obtained 
with the axial HFB code has a prolate deformation with 
/3 = 0.143, and both neutrons and protons are supercon¬ 
ducting with A„ = 0.702 MeV and Ap = 0.517 MeV. The 
experimental averaged pairing gaps are = 

1.170 MeV and A^^^^oxe) = 1.014 MeV. 

The response functions computed from the FAM are 
S{Nr,,Nn) = 12.704 MeV-\ SiNr^Np) = 8.725 MeV^, 
and S{Np,Np) = 3.083 MeV“^. We took a numer¬ 
ical average of S{Nn,Np) and S{Np,Nn) for the off- 
diagonal term. The Thouless-Valatin moments of inertia 
are Jnn = 11.986 MeV“^, 57„p = —4.236 MeV“^, and 
Jpp = 2.909 MeV“^. The opposite sign of Jnp is consis¬ 
tent with the isorotation picture, whose rotational energy 
is proportional to T(r -|- 1) [7, 65, 66], and produces a 
negative sign for the neutron-proton term. 

Figure 6 shows the pairing rotational energies mea¬ 
sured from ^^°Xe along the Xe isotope direction, the 
N = 76 isotone direction, the A = 130 isobar direction, 
and the = 11 direction. Clear parabola patterns are 
seen in the pairing rotational energy from the HFB and 
the experimental data in Fig. 6 a) and b), although the 
reference state has a prolate deformation and the intrinsic 
shape changes nucleus by nucleus. The Thouless-Valatin 
inertia explains the neutron and proton pairing rotational 
energy in the vicinity of the reference state. 

In Fig. 6 c), the pairing rotational energy along the 
A = 130 isobar is shown. Again a parabola pattern 
is found along the isobar, and both the HFB and the 
Thouless-Valatin inertia explain the experimental data. 
From Eq. (74) the pairing rotational energy along the 
isobar is given by 


Erot{N, Z) 


1 


2Jnn 2X, 


2Jpp) 


(79) 


with ATz = 11 —T^. This represents the isorotational en¬ 
ergy which restores the broken isospin symmetry. The co¬ 
efficient in front of (AT^)^ is 0.44 MeV, and explains the 
systematic behavior of the binding energies. For compar¬ 
ison we also show the pairing rotational energy with the 
Thouless-Valatin inertia but without the neutron-proton 
term in the same figure. The value is close to the Belyaev 


inertia that does not contain the neutron-proton term 
either, indicating the importance of the neutron-proton 
term for the isorotation. 

The last example is the = 11 nuclei shown in Fig. 6 
d). All the calculations and the experimental data show 
the pairing rotational energy is small. This degree of 
freedom is associated with the breaking of global gauge 
invariance (the total particle number symmetry) [65]. In 
fact, the pairing rotational energy from the Thouless- 
Valatin inertia along a constant Tz line is written as 


Erot{N,Z) 


2Jn 




(80) 


and the coefficient in front of (AA)^ = (A—130)^ is — 0.02 
MeV. Our functional preserves an approximate isospin 
symmetry well, and only the Coulomb term breaks it. 
Using independent pairing strengths for neutrons and 
protons introduces explicit isospin symmetry breaking 
in the pairing channel, and may generate pairing rota¬ 
tional energy associated with global gauge invariance. 
The quadrupole shape changes rapidly along the = 11 
chain, and the correlation energy from the deformation 
is larger than the pairing rotational energy. This causes 
a deviation from the parabola curve in the pairing rota¬ 
tional energy. 

As a whole, the pairing rotational description based on 
the open shell nucleus ^^°Xe works well, and explains the 
experimental binding energy systematics of neighboring 
even-even nuclei around ^^°Xe. This indicates that the 
intrinsic superconducting HFB state of ^^°Xe contains 
information of neighboring even-even nuclei. With this 
generalization of the pairing rotation to the neutron and 
proton mixed modes, we can take an arbitrary supercon¬ 
ducting nucleus as a reference state of the pairing rota¬ 
tion. Systematic analysis of the pairing rotational modes 
and moments of inertia in this direction is in progress 
[67]. 


VIII. CONCLUSION 

We formulated linear response theory in the presence 
of zero-energy NG modes. We showed that the Thouless- 
Valatin inertia of the NG mode is derived from the zero- 
frequency response to the momentum operator of the NG 
mode. Combined with the finite-amplitude method for 
nuclear energy density functional theory, we can compute 
the Thouless-Valatin inertia very precisely and efficiently. 
This formulation also provides the expression for the con¬ 
jugate coordinate operator that will be necessary for the 
spurious mode removal procedure of the FAM. The for¬ 
mulation was numerically tested in detail for the cases of 
the trivial center of mass mode and the pairing rotational 
mode. Although the center of mass mode is not at zero 
energy in practical numerical calculations, the strength 
function at zero frequency provides a precise Thouless- 
Valatin inertia. The realistic applications for the neutron 
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FIG. 6. (Color online) a) Nentron pairing rotational energy along Xe isotopes, b) proton pairing rotational energy along N = 76 
isotones, c) neutron and proton pairing rotational energies along the A = 130 isobar, and d) along = 11 nuclei measured 
from ^®°Xe (N = 76, Z — 54). The experimental neutron and proton chemical potentials are An,exp = —7.926 MeV and 
Ap.exp = —7.149 MeV. The red solid curves are the pairing rotational energies with the Thouless-Valatin inertias in Eq. (74), 
black dotted curves are the same energies but without the neutron-proton terms. The blue dashed curves are the pairing 
rotational energies from the Belyaev inertias without neutron-proton terms. 


and proton pairing rotational moments of inertia are pre¬ 
sented for Sn isotopes and N = 82 isotones, respectively. 
We then consider the situation where both neutrons and 
protons are in a superconducting phase. We show that 
the neutron and proton pairing rotational degrees of free¬ 
dom are mixed in the QRPA normal inodes, and the 
neutron-proton term in the pairing rotational energy is 
generated. As a realistic application, we show that the 
ground state energies around ^^°Xe can be interpreted in 
terms of the pairing rotation picture. 

One interesting future application of this formalism for 
NG modes is in the computation of the Thouless-Valatin 
rotational moment of inertia. A systematic comparison 
of the energy of the 2f states of deformed nuclei with 
experimental data may clarify the property of the un¬ 
constrained time-odd term of the nuclear EDF. Once im¬ 
plementation into the symmetry unrestricted HFB code 
has been completed, this approach can serve as an ef¬ 
ficient technique for deriving the Thouless-Valatin mo¬ 
ments of inertia of three-dimensional collective rotation 
for the five-dimensional quadrupole collective Hamilto¬ 
nian [16]. 


Deriving the pairing collective Hamiltonian [68-70] 
based on the nuclear EDF and discussing the anhar- 
monic and large-amplitude aspect of the collective pair¬ 
ing motion and the coupling of the pairing vibration and 
pairing rotation will be challenges for the future. Ex¬ 
tension of the neutron and proton pairing rotations to 
full SU{2) isorotation by including the neutron-proton 
pairing within the isospin-invariant nuclear density func¬ 
tional theory [71-73] is another future challenge for un¬ 
derstanding the role of the isospin symmetry in low- 
energy nuclear collective motion. 
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